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ABSTRACT 

A grid of numerical simulations of double-diffusive convection is presented for astrophysical conditions. As in laboratory 
and geophysical cases convection takes place in a layered form. A translation between the astrophysical fluid mechanics 
and incompressible (Boussinesq) approximation is given, valid for thin layers. Its validity is checked by comparison 
of the results of fully compressible and Boussinesq simulations of semiconvection. A fitting formula is given for the 
superadiabatic gradient as a function of this parameter. The superadiabaticity depends on the thickness d of the double 
diffusive layers, for which no good theory is available, but the effective He-diffusion coefficient is nearly independent of 
d. For a fiducial main sequence model (15 Mq) the inferred mixing time scale is of the order 10^° yr. 
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1. Introduction 

In models of stellar structure, situations are found where 
the heavier products of nuclear burning provide stability 
to a zone which otherwise would be unstable to convective 
overturning. Such a zone, or part of it, would become con- 
vective if something managed to mix its composition (R.J. 
Tayler, 1953). The question whether such a zone should be 
treated as if it were mixed or not has become known as 
the semiconvection problem. Answers to this question dif- 
fer substantially. In practice, recipes are used containing a 
free parameter that allows the degree of mixing to be var- 
ied. Calculations in which such a parameter is adjusted to 
match observations are then called 'with semiconvection'. 
Commonly used prescriptions are those of Langer (1985) 
and Maeder (1997). 

The presence of a semiconvective zone has only a minor 
effect on the thermal structure of the star. The assumed 
amount of mixing of composition is critical, however, be- 
cause the evolution of the star is sensitive to the precise 
distribution of products of nuclear burning with depth in 
the star. The main goal of a theory for semiconvection is 
thus a good determination of the rate of mixing. From the 
perspective of the stellar evolutionist, the theory would ide- 
ally provide a formula for the rates of mixing and energy 
transport (the effective diffusivities), as functions of local 
thermodynamic state and composition, and their gradients. 

In Spruit (1992, hereafter S92) such formulas were de- 
rived, adapting the known physics of double-diffusive con- 
vection to the case of a stellar interior. In the following, 
numerical simulations are used to measure mixing rates 
and their dependence on astrophysical conditions. They are 
compared with S92, and used as a basis for updated fitting 
formulas for the mixing rates in semiconvective zones of 
stars. 
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2. Semiconvection and double-diffusive convection 

Situations where a fluid is stabilized by the density gradi- 
ent due to a dissolved heavy constituent occur in nature. 
An example is convection under the arctic ice sheet (cooled 
from above, stabilized by the salts dissolved in the sea wa- 
ter). Particularly intensively studied are East- African rift 
lakes (lakes Kivu, Nyos and Mounon, cf. Schmid et al. 
2010). These are heated from below by volcanic activity, 
which also is a source of dissolved gases (carbon dioxide and 
methane, hereafter the 'solute'). Their density stratification 
is stabilized against convection by the stable gradient re- 
sulting from the weight of the carbon dioxide. Efforts to pre- 
vent catastrophic release of carbon dioxide (Lake Nyos, e.g. 
Sigvaldason 1989) or commercial exploitation of methane 
(Lake Kivu, Nayar 2009) have led to extensive study of 
the fluid flows, heat flux and mixing rates in these natural 
double-diffusive systems. 

The gradients in temperature and solute in these lakes 
are observed to be 'stepped': consisting of a stack of thin 
layers (decimeters to decameters). Inside a layer, overturn- 
ing convection keeps the composition nearly uniform, with 
stable jumps in temperature and composition separating 
the layers. The physics involved is easily reproduced under 
controlled laboratory conditions (Turner 1985) The lay- 
ers are very long-lived: of the order of months or more in 
the geophysical examples mentioned, orders of magnitude 
longer than the convective turnover times inside the layers. 

In the stable steps between the layers the transport of 
the stabilizing solute takes place by diffusion instead of con- 



^ Also on a coffee table. A latte macchiato in a tall glass often 
shows the effect nicely. After the coffee is added to the milk, a 
stably stratified gradient of milk/coffee mix develops (showing 
internal gravity waves in the form of a sloshing motion with 
a period of a few seconds). After about a minute, the initially 
smooth gradient starts dividing into thin (a few mm) layers, 
visible at low contrast. In the course of several minutes these 
merge into a smaller number of more clearly defined layers. 
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vection. This strongly limits the effective transport of so- 
lute through the double-diffusive stack of layers. Residence 
times on the of order 1000 yrs are inferred for the solutes in 
lake Kivu, for example (Schmid et al. 2010). This is 8 orders 
of magnitude longer than the convective turnover times in 
these layers. The transport of heat is also strongly reduced, 
this is exploited for heat storage in solar ponds (cf. Lu and 
Swift, 2001). 

Theoretically, the observed layered nature of double- 
diffusive systems is well understood. Central to this un- 
derstanding is the fact that linear stability analysis does 
not provide relevant clues to their behavior, because the 
double-diffusive case of thermal convection stabilized by a 
slowly-diffusing solute is subcritical. That is, stable over- 
turning flows occur at a temperature gradient below the 
onset of linear instability. 

Linear instability predicts internal gravity modes to set 
in above some critical value of the temperature gradient, 
growing in amplitude by the effect of thermal diffusion, the 
so-called Kato oscillations (Kato 1966). Such oscillations 
(cf. movie at Fig. pj transport a negligible amount of heat 
or solute, compared with overturning motions of the same 
amplitude. For this reason alone, linear stability arguments 
cannot be used for useful estimates of the mixing rate in 
semiconvective zones. More important is the subcritical na- 
ture of double-diffusive convection. Proctor (1981) shows 
that, in the limit of vanishing diffusivity of the solute, the 
layered form of convection can occur whenever the temper- 
ature gradient is unstable to convection (i.e. 'according to 
Schwarzschild'), irrespective of the strength of the stabiliz- 
ing component. 

The reason for this subcritical behavior can be under- 
stood with an energy consideration. The amount of energy 
it takes to overturn a layer of thickness d against a stable 
gradient scales as (P. Per unit mass the expense in initial 
energy needed to put the system into its finite-amplitude, 
layered state thus vanishes as d, down to the value where 
energy loss by viscous damping stabilizes the system. This 
minimum is set by the critical Rayleigh number for ordinary 
convection in a layer of thickness d. A small initial pertur- 
bation, or perhaps an initial Kato oscillation, is sufficient 
to provide the energy for overturning into thin layers. Once 
established, this layered state is a stable form of convection. 
This agrees with the observation that in laboratory exper- 
iments and geophysical systems like lake Kivu mentioned 
above the layering first sets at a small thickness (cf. the 
latte macchiato experiment in the footnote above). 

2.1. Semiconvection 

It is sometimes argued that the geophysical and laboratory 
examples of double diffusive convection cannot be applied 
to a stellar interior because of different physics. At the level 
of physics in evidence in current recipes for semiconvection, 
however, there is no distinction between these cases. The 
equations of fluid dynamics used are the same. In contrast 
with astrophysics, an incompressible approximation is often 
used in geophysics, but this is not mandatory. The com- 
pressibility of water can be, and has been included widely 
in geophysical fluid mechanics. 

The fluid is described by the thermodynamic variables 
defining its local state (e.g. pressure and density), the gra- 
dients of these with depth, the transport properties that are 
determined by the thermodynamic state (viscosity v, ther- 



mal diffusivity kt and solute diffusivity Kg), the accelera- 
tion of gravity, and the equation of state. Taken together, 
these quantities form a large parameter space, and it might 
be concluded that realistic numerical simulations of semi- 
convection would have to be done for individual zones in 
individual stellar models. 

The equations of fluid dynamics have symmetries, how- 
ever, so that the independent degrees of freedom are far 
fewer. They can be represented by 5 dimensionless parame- 
ters: a Rayleigh number Ra, the layering thickness e = d/H 
in units of the pressure scale height Hp, the Prandtl number 
Pr= i^/kt, the Lewis number Le = kq/kt, and a density 
ratio Rp which measures the ratio of stabilizing (solute) to 
destabilizing thermal gradient. The behavior of semicon- 
vection at any point in a star can be deflned in terms of 
these parameters. 

This is discussed in more detail in section [3j where it 
is also shown how in the limit e ^ 1 the astrophysical 
problem can be translated into an equivalent incompressible 
(Boussinesq) problem. In section 5.4 we verify this with a 
direct comparison between results from fully compressible 
and Boussinesq simulations. In this limit, the e disappears 
from the problem, reducing the number of parameters from 
5 to 4. 

By a fortunate coincidence, it turns out that the results 
are effectively independent of Pr, as long as it is small. 
Which is in fact the case in a stellar interior. This further 
reduces the number of independent parameters to only 3. 
Since measurement of the mixing rate in each individual 
case does not require a very expensive simulation, a sig- 
nificant volume of astrophysically relevant parameter space 
can be covered. 

The semiconvective zone in a stellar model is sufficiently 
limited in extent and consequences that the overall struc- 
ture of the star does not depend much on the way semicon- 
vection is calculated. The additional diffusion of Helium 
by semiconvective mixing can be important for later evolu- 
tionary stages, but during the semiconvective phase itself 
the thermal structure of the star is not affected much. The 
consequence of this is that in contrast with to laboratory 
and geophysical situations, in a stellar model the heat flux 
F can be considered as known, rather than the temper- 
ature gradient. Since the radiative contribution F,- to the 
heat flux is known to good approximation from the thermal 
structure of the star, the convective heat flux Fc = F — F^ 
transported by semiconvection is then also known. What is 
not yet known is the efficiency of convection: i.e. how close 
to the adiabatic gradient the mean thermal gradient will 
turn out to be as a result of semiconvective transport. 

2.2. Layered convection 

Asymptotically at high Rayleigh number, the convection 
inside each layer of the semiconvective stack is driven en- 
tirely by boundary layers at the top and bottom steps of the 
layer, much like in laboratory convection in a box. Except 
for these thin boundary layers, the (horizontal averages of) 
entropy and composition are almost uniform inside the lay- 
ers. 

Under the conditions in a stellar interior, the thermal 
diffusivity is much larger than the diffusivity of the so- 
lute (Helium in Hydrogen, say). The thickness of the so- 
lute boundary layers is then much smaller than the thermal 
boundary layers, and the convective flow almost the same 
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as in the absence of a stabilizing solute. Under these asymp- 
totic conditions, the dependence of heat flux on Rayleigh 
number can be taken from a simple estimate, as was done 
in S92, or a laboratory result can be used (e.g. Niemela et 
al. 2000). With the convective flow thus known, the flux of 
solute can be calculated as well. 

One of the consequences of such an asymptotic model 
is the prediction that the effective diffusivity of the solute 
scales with the square root of its microscopic value. This is 
the well-known relation derived in various ways for double- 
diffusive convection (cf. Turner 1985), and verified in labo- 
ratory experiments. It will show up again in the numerical 
experiments reported below. 

2.2.1. Layer thickness 

The main uncertainty in applying results on double- 
diffusive convection to stars (or any other double-diffusive 
situation) is the poorly known physics which determines 
the thickness of the individual layers in the stack. The layer 
thickness thus remains a free parameter of the problem. 



2.3. Boundary layers 

Each of the layers in the stack is separated from the next 
by a step in composition and temperature. The step con- 
sists of three nested boundary layers, for temperature, com- 
position, and flow speed, called the thermal, solute and 
viscous boundary layers respectively. Their thicknesses are 
governed by the thermal diffusivity kt, solute diffusivity Kg, 
and (kinematic) viscosity i'. The highest diffusivity (ther- 
mal) has the widest boundary layer. For astrophysical con- 
ditions, kt ^ Ks, the solute and viscous sublayers are so 
thin compared with the thermal boundary layer that their 
influence on the convective flow pattern is negligible. The 
viscous sublayer is different from the other two. Whereas 
temperature and solute change rapidly across the step, the 
parallel flow component varies smoothly across the step. 
The stable step in composition conflnes the flow to the 
layer. Apart from of deformations of the step in the form 
of stable interface waves, the vertical velocity component 
vanishes at the step. 

Summarizing the above, in the astrophysical limit kt 3> 
Kg, I/, the flow in the layer behaves to a good approxima- 
tion like convection between horizontal plates, with free slip 
conditions at these boundaries. This allows us to make esti- 
mates of the thickness of the boundary layers. The bound- 
ary layers determine the fluxes of heat and solute through 
the layer, i.e. the numbers we are interested in. The need to 
properly resolve them numerically determines the required 
number of grid points in the vertical direction. 

Let Tc be the time the overturning flow is in contact 
with the boundary, before descending/ascending into the 
interior of the layer. Diffusion of heat through the step over 
this time sets the thickness St of the thermal boundary 
layer: 



(ktTc) 



1/2 



(1) 



Assuming that the convective cells driven by the boundary 
layers have a horizontal width of the order of the thickness d 
of the layer, Tc is of the order of the convective overturning 
time, Tc ~ d/vc, where the convective velocity Vc can be 



estimated by the mixing-length expression, 

d^ 



V 



3(V-Va) 



H' 



(2) 



[Standard notation: pressure scale height Hp — 
l/(dlnP/dz), V = dlnT/dlnP, Va the corresponding adi- 
abatic gradient, and rj = — (91np/91nT|p is a number of 
order unity (called 6 in Kippenhahn & Weigert 1990, p39), 
?7 = 1 in a fully ionized ideal gas] . Define a modified Rayleigh 
number Ra* by 

Ra*=PrRa. (3) 
In astrophysical notation, it is given by 



Ra* - |:(V-Va)/(/3)^ 

11 /S^rp 



,2^2 



(4) 



where g is the acceleration of gravity, and /(/3) a factor 
depending on the ratio, /3 — Pg/P, of gas pressure Pg to 
total (gas plus radiation) pressure, P = Pg + Prad- Ra* 
thus contains only the thermal diffusivity, not the viscos- 
ity. Except for Rayleigh numbers close to the critical value 
Rac for instability to convection, it is the quantity that de- 
termines the heat fiux (Nusselt number), in the limit of low 
Prandtl number. The simulations below will also address 
conditions close to marginal, but Ra* is the meaningful 
number for translating them to astrophysical application, 
where Ra ^ Rac. Using equations ([2| and ([l]). 



4 



r] Ra* 



(5) 



Since 8f /rj is a numerical factor of order unity, our estimate 
of the thermal boundary layer thickness is 



6t /d ~ Ra^ 



-1/4 



(6) 



The Nusselt number, the ratio of fiux in the presence 
of the flow to the fiux in the absence of a flow is then 
just the ratio of the gradient in the boundary layer to the 
average gradient between top and bottom. The expected 
dependence on Ra* is thus, with (|6]): 



Nut ~ d/S^ ~ Ra, 



1/4 



(7) 



The thickness of the solute (Helium-) boundary layer 
can be estimated by the same reasoning. 



1/2 



Lei/^(KTT,)i/2=Lei/2^T, 



(8) 



where Le is the Lewis number Le = kq/kt as before. Being 
the thinnest boundary layer in the problem, Sg determines 
the numerical resolution needed. Since the fine structure in 
the interior of the layer largely consist of plumes 'peeled 
off' from the boundaries, the same resolution is needed in 
the interior of the layer as well, and there is no need or 
justification for using non-uniform grids. 

2.4. Density ratio, mixing rate 

The problem also depends on the relative strength of the 
stabilizing solute gradient relative to the destabilizing ther- 
mal gradient. This can be measured in terms of the thermal 
and solute buoyancy frequencies iVr, Ng: 



n4 



H 



(Va-V), 



(9) 
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(10) where g is a constant depending on the local thcrmody- 



where /i is the 'mean molecular weight'. The density ratio 



Ro 



(11) 



is then the dimensionless measure we will use for the rela- 
tive strength of the stabilizing solute gradient. The minus 
sign is added for convenience to the definition, since < 
for a convectively unstable stratification. Semiconvection, 
i.e. a stratification 'between Schwarzschild and Ledoux' 
then corresponds to i?p > 1. 

The amount of solute transported across the layer is 
limited by the fact that convective plumes develop only 
from material with a net buoyancy of the unstable sign. 
This limits the density contrast of the stable solute carried 
to a fraction 1/Rp of the density contrast across the layer. 
With the same reasoning used above for the heat flux, this 
gives a relation between the solute Nusselt number Nus and 
Nut : 

Nus = Le~^/2NuT/i?p. (12) 

(cf. S92). Solute with a higher density contrast is not mixed 
into the convective flow: it forms a stagnant region around 
the interface. This region has a width 



SsR, 



(13) 

which decreases the solute gradient and thereby the solute 
diffusion rate at the interface. The effect is equivalent to 
the reduction of Nus by the factor Rp in (12). 

In the above (as in S92), we have implicitly assumed 
that Rp, while > 1, is still low enough that the stagnant 
region is thinner than the thermal boundary layer thickness: 



ds < St, 



(14) 



so that it does not interfere with the convective flow. An 
improvement for the more general case yields 



Nus = (1 + Le"i/V^p)NuT 



(15) 



As discussed in section [673| the Lewis number is so small in 
a stellar interior that this correction is not likely to become 
important, and estimate ( 12 ) is still good enough. The heat 



flux can be considere d as flxed by the stellar structure, so 
Nut is fixed and Eq. ( 12 1 also fixes the mixing rate, largely 
independent of the unknown layer thickness d. 



3. Boussinesq limit: thin layers 

In a stellar interior, radiation carries a heat fiux even in a 
convectively neutral stratification. This requires some care 
when making the connection with the Boussinesq case. 



3.1. heat flux 

Let F^, be the radiative and hydrodynamic (semicon- 
vective) contributions to the heat flux in the star, where the 
superscript designates quantities under the compressible- 
fluid conditions in a star, while ^ will be used for the cor- 
responding quantities in the Boussinesq (incompressible) 
model. The radiative heat flux is proportional to the tem- 
perature gradient V: 



namic state, and F^^, F^^ are the contributions to the 
radiative heat flux of the adiabatic and the superadiabatic 
parts of the mean temperature gradient. In the Boussinesq 
model, the contribution i^ra is absent: the convective and 
radiative heat fluxes F^ , F^ are governed by the same tem- 
perature gradient. Related to this, the Boussinesq model 
has one fewer parameter: the pressure scale height H. 
Because of this difference, the heat flux in the stellar (com- 
pressible) model cannot be compared directly with the heat 
flux in an incompressible model. Instead, the ratio of 
convective to radiative heat flux in the Boussinesq model is 
to be identified with the ratio P — F^/ F^^ of the convective 
fiux to the superadiabatic component of the radiative flux 
in the star, in the limit H — >■ oo. If e = d/H, where d is the 
thickness of the convecting layer, then for a given value of 
the modified Rayleigh number Ra* (cf. Massaguer & Zahn 
1980): 



/B(Ra,;M) = f (Ra„eiO;M), 



(17) 



where {p} stands for the other dimensionless parameters of 
the problem: Pr, Le and Rp. The Boussinesq model thus is 
the limit e 4- taken at fixed Ra* . If the semiconvective layer 
thickness e is 'sufficiently small' (in some sense be verified 
by numerical tests) , the compressible case can be comp are d 
directly with the Boussinesq model, in the sense of eq. (17). 

As mentioned in the introduction, this makes the semi- 
convection problem particularly amenable to numerical 
simulation. In contrast with a convective stellar envelope, 
for example, with its many scale heights to be covered, the 
layered nature of double diffusive convection puts it in a pa- 
rameter range that is much more accessible with realistic 
ab initio calculations. 

To complete the Boussinesq case, we s till n eed a trans- 
lation of the density ratio Rp. With eq. (Ill]) this can be 
done uniquely in terms of the buoyancy frequencies, see 
eqs. ( [22p3l ) below. 



4. Numerical simulations 

All equations were calculated on a 2D rectilinear Cartesian 
grid in terms of finite differences. The set of equations 
are implemented into the ANTARES software framework 
(Muthsam et al. 2010). Advective currents are solved by 
a weighted essentially non-oscillatory finite volume scheme 
in fifth order (Shu & Osher 1988), the physical diffusion 
is handled by a fourth-order finite difference discretization. 
A second order total variation diminishing scheme as time 
integrator is chosen. To avoid odd-even decoupling, a MAC 
grid, which locates vector variables at cell faces and scalar 
variables at cell centres, is used. 

For a detailed description of the numerical solution of 
the binary mixture equations presented here see Zaussinger 
(2011). 

4.1. Boussinesq approximation for a binary mixture 

In the Boussinesq approximation the continuity equation is 
reduced to that of an incompressible fluid. 



qV = gVa + <z(V - Va) = F:^ 



(16) 



dp 
di 



(18) 
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such that in the equation of motion the density is taken to 
be a constant poi 



du 

Id 



1 

Po 



VP + {-aTe + asS)g + V -{uVu), (19) 



where the 'expansion coefficients' ax, as describe the den- 
sity effects of variations 8, S in temperature and solute, as- 
sumed small compared with the mean temperature O and 
solute S. They are given by advection-diffusion equations: 



de 

df 

d5 
'dt 



(20) 



(21) 



where v is the kinematic viscosity. The mean gradients V8, 
VS" can be expressed more usefully in terms of the buoyancy 
frequencies: 

7V| = ax g • Ve, (22) 



Nl 



as g • VS, 



(23) 



The density ratio Rp is then defined as in the compressible 
case, eq. dTTl). 



The numerical algorithm solving this set of equations is 
based on a semi-implicit scheme. Intermediate values for the 
velocity field u* are calculated explicitly from the equations 
of motion. By the nature of the incompressible equations, 
the pressure update is done implicitly, by solving a Poisson 
equation: 



AV ' 



(24) 



The resulting pressure P leads to the required diver- 
gence free velocity field at the new time step n -|- 1. 



=u* - —VP. 
Po 



4.2. Compressible fluid equations for a binary mixture 



(25) 



Verification of the Boussinesq results has been done with 
simulations of the fully explicit compressible fluid equa- 
tions. The fluid is assumed to be an ideal gas, which is a 
good approximation to a binary gas mixture of our interest. 
These are the continuity equation 



dp 
dt 



V • (pu) = 
the partial density equation 

V • {pcu) — V • {pKc^c) 
the momentum equation 

V • {puu + PI) = pg^ + V • T 



d{pc) 
dt 



d{pu) 
dt 



the total energy equation 

de 
dt 

and the equation of state 



V • [u{e + P)] = p{gu) + V • (wr) 



(26) 



(27) 



(28) 



(29) 



P = 



TZpT 

p{c) 



(30) 



where c is the solute mass fraction, Kc is the solute diffu- 
sion coefficient, r is the viscous stress tensor and p — p{c) 
the mean molecular weight. 



4.3. Units, boundary and initial conditions 

As unit of length we use, for the Boussinesq cases the layer 
thickness d, for the compressible calculations the pressure 
scale height H. The nominal convective turnover time d/vc 
from ([2J) is used as unit of time. As units of temperature 



(potential temperature Q in the compressible case) we use 
1/aT, for solute concentration 1/as. The density ratio then 
becomes Rp = VS/VT. The Boussinesq equations are in- 
variant to arbitrary additive constants in temperature and 
solute. We set these such that T and S are zero at the top 
boundary. 

Because of the symmetries of the problem, there are 
fewer independent parameters than physical variables de- 
scribing it. Hence some of the physical quantities appearing 
in the problem can be set to unity. We choose for these: 
the temperature difference between top and bottom of the 
layer, the density at the bottom of the layer and the accel- 
eration of gravity. The Rayleigh number Ra* is then con- 
trolled through the thermal diffusivity kt, the solute differ- 
ence between top and bottom though the density ratio Rp, 
and the solute diffusivity through the Lewis number Le. 

Most of the calculations were done in a box simulating a 
single layer from the double-diffusive stack, so the top and 
bottom boundaries coincide with the steps between lay- 
ers. This ignores the distortions of the interfaces by surface 
waves, but since the essence of the double layering phe- 
nomenon is that the transport across the interface is by 
diffusion, this is not expected to make a big difference. To 
check that this is indeed the smaller set of simula- 

tions wa s do ne in which a step is present inside the volume 
(section 5.5 1. The vertical boundary conditions are thus 
taken to be impermeable and stress-free. In the horizontal 
direction periodic conditions are used. 



Uz 


= 


z = 


0,1 


(31) 


dux 








(32) 


dz 





z = 


0,1 


S 


= Rp 


z — 





(33) 


T 


= 1 


z = 





(34) 


S 


= 


z = 


1 


(35) 


T 


= 


z = 


1 


(36) 



As initial condition the stratification of temperature and 
solute was taken to be horizontally uniform with either 
a linear gradient between the values at top and bottom 
(the 'linear' case below), or something approximating the 
boundary layer structure expected of the final state (the 
'step' case). Small initial random perturbations are applied 
on the solute field. 

The numerical algorithm for setting up the initial con- 
ditions for the compressible fluid equations is an extension 
of the procedure presented in (Muthsam et al. 1995, 1999). 
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Pr=10"' and Le=10"^ 



Fig. 1. Flow structure in a double diffusive layer. The 
temperature field (top) is more diffuse than the solute 
('Helium') concentration, as a result of the high thermal 
diffusivity . 



4.4. Numerical setup 

Being the thinnest boundary layer in the problem, deter- 
mines the numerical resolution needed near the boundaries. 
Since the fine structure in the interior of the layer largely 
consist of boundary layers 'peeled off' from the boundaries, 
the same resolution is needed in the interior of the layer 
as well, and there is no need or justification for using non- 
uniform grids. 

Using Ra* = 1.6 x 10^ and a Lewis number of Le = 0.1 
would lead to solute boundary thickness of 1.6% of the layer 
thickness d. With the high-order spatial discretization used, 
a minimum resolution of 3 points is then needed in the 
boundary layers, which translates to 200 points in vertical 
direction for this case. Convergence tests showed that a nu- 
merical resolution of 300 points was sufficient for the range 
in Lewis number considered. Most of the calculations were 
done with a horizontal-to- vertical aspect ratio of 2:1. A few 
tests with different ratios showed that this choice does not 
affe ct the measured Nusselt numbers significantly (section 
5^ 



5. Results of numerical simulations 

The numerical results presented here are based on about 
100 numerical simulations done in the Boussinesq approxi- 
mation and about 20 simulations performed with the fully 
compressible code. 

An example is shown in Fig. [l] with Pr — 0.1, Ra* = 
5x 10^, Le — 0.01, Rp=1.15. It shows the key characteristic 
of double diffusive convection: the boundary layers and the 
plumes of the solute are narrower than those of tempera- 
ture, on account of the low solute diffusivity. The fiux of 
solute carried by the velocity field is correspondingly lower, 
and vanishes in the limit Le — (cf. eq. 38). 



R =1.15 
R =1.50 
R =2.00 
Rj=5.00 
Ra.=16000 
Ra.=50000 
Ra.=1 60000 
Ra.=500000 
Ra.=1 600000 




NUj-1 



Fig. 2. The relation between Nut and Nus for Pr = 10^^ 
and Le = 10~^. Parameter along the solid curves is the 
Rayleigh number Ra», colors indicate density ratio Rp. 



Dotted lines show the theoretical estimates from eq. 38 



For comparison with the numerical results, we use the 
theoretical estimates in |2.3| . For large Nut, Nus, the anal- 
ysis in 2.3 predicts a classical double-diffusive square root 



dependence on the ratio of diffusivitics: 



Nus « (1 + Lc"i/V^p) Nut 



(37) 



which is observed in laboratory experiments. This estimate 
assumes the boundary layers to be thin compared with the 
layer thickness. An estimate which is slightly better at low 
Nusselt numbers is obtained by noting that the convective 
flux F~Fo is a more relevant measure of the double diffusive 
transport efficiency than the total flux F; this yields: 



(Nus - 1) « (1 + Le-'/^/Rp) (Nut - 1). 



(38) 



We use these estimates for comparison with the numerical 
results. 



5.1. Dependence on Ra^ and Rp 

Figs.[2]shows the dependence of the measured Nusselt num- 
bers Nut and Nus on the parameters Ra* and R^, for the 
case Le= 0.01, Pr=0.1. The results show considerable scat- 
ter, as expected from the limited number of overturning 
times for which the simulations were run. Nevertheless, they 
appear consistent within about a factor of 2 with theoretical 
estimate (38). 



The numerical results thus conflrm that the ratio of 
thermal to solute transport does not depend much on pa- 
rameters other than Lewis number and density ratio. The 
actual transport efflciency of both of course does depend 
on the Rayleigh number Ra* . This is shown in Fig. [3j and 
compared with the thermal convection results discussed 
below in section 16.21 In S92 an estimate of this relation 
was made for the asymptotic dependence at large Ra*: 
Nut ~ 0.5 Ra*^/'' (for the essence of the derivation see 
section 2.3). As above, the estimate can be made more ac- 
curate at lower Nusselt numbers by subtracting the diffusive 
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Fig. 3. Nut as function of Ra* for Pr — 10 
10-2. 



and Le = 



flux: 



1 w 0.5 Ra,^/^. 



(39) 



It is seen that this fits the numerical Nusselt numbers 
for density ratios near unity, but overestimates the heat flux 
at higher Rp. For the relation between Nus and Nut this 
has little effect, since both are similarly affected by Rp. As 
discussed below, since the heat flux is fixed by the structure 
of the star, this uncertainty should not affect the expected 
mixing rate in a stellar interior much. 

The same data is shown again in Fig. |4j with Rp on the 
horizontal axis and Ra* as parameter. 
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Fig. 4. Nut as function of Rp for Pr = 10 ^ and Le = 10 



5.2. Dependence on aspect ratio 

The aspect ratio of 2 : 1 used in the results reported above 
was tested against 5 : 1 and 10 : 1 at the same spatial res- 
olution. For the reference simulation (Pr = 0.1, Le = 0.01, 
Ra* = 10^ Rp = 1.15) we find Nug = 90 and Nut = 8.5. 




Fig. 5. Example of the development of an overturning 
flow from Kato oscillations. Time f rom left to right 
and top to bottom. See also movi e at http://www.mpa- 
|gar ching . mpg .de/^henk/ movie . avi 



By comparison the simulation with 5 : 1 results in Nus — 90 
and Nut = 9.0. The most extended box with an aspect ra- 
tio of 10 : 1 and a spatial resolution of 1500 x 300 has 
Nusselt numbers of Nus = 80 and Nut = 8.75. The aspect 
ratio thus has no signiflcant influence on the dependences 
of the fluxes on input parameters in our simulations, within 
the fluctuations due to the stochastic nature of the flow. 



5.3. Dependence on initial stratification 

As initial state we used either a constant linear gradient 
of temperature and solute between top and bottom val- 
ues ('linear'), or a profile with boundary layers of width 
as estimated in 2.3 ('step'). The linear case shows how the 
oscillatory phase due to th e K ato instability develops into 
overturning flow, see Figs. |5|6[ The duration of the initial 
formation process is mainly determined by Rp and Ra*. 
The end state in both cases is statistically the same. The 
'step' as initial condition also covers cases that are stable 
in linear theory because of the subcriticality of the system, 
and can gain a large factor in computing time for small 
values in Ra* and high Rp. It was used in all the results 
reported above. 



5.4. Comparison witli compressible results 

The compressible simulations are based on a 5th order 
weighted ENO scheme. Compressible fluids lead to restric- 
tions in time stepping (due to the need to resolve sound 
waves). The compressible simulations take up to 100 times 
longer for the same resolution compared to the incompress- 
ible solver. Therefore only a few tests have been done with 
the compressible code. The degree of compressibility is gov- 
erned by the e = d/H of layer thickness to pressure scale 
height; in the limit e — >■ there is a direct translation be- 
tween the compressible and the Boussinesq case (section 

t. The results of a numerical comparison with e = 0.1 is 
own in Table 1. The resulting Nusselt numbers do not 
differ significantly. Simulations done with e = 1.0 behave 
quite similar to these done with e = 0.1. The mixing pro- 
cesses do not significantly differ as long as the Rayleigh 
number is high enough, Ra > 5.0 x 10^. At the present 
level of accuracy (a factor of 2, say), we conclude that the 
Boussinesq approximation gives the right results even for 
layer thicknesses approaching a scale height. 
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Fig. 6. Top panel: evolution of the Nusselt number for the 
linear initial stratification. Convective cells get established 
from Kato oscillations (cf. Fig. [5]) after about 10 turnover 
times. Starting the simulation irom a step (bottom) saves 
computing time. At the end of the runs the Nusselt num- 
bers of both simulations are the same within the statistical 
variations. 
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33 


11 


26 


10 



Table 1. Comparison of compressible and incompressible 
simulations. Thermal (t) and solute (s) Nusselt numbers 
from the Boussinesq (^) and compressible {^) results. Layer 
thickness d is 0.1 pressure scale height. 



5.5. Multi-layer simulations 

In all of the above wc have assumed that the interfaces 
between the double diffusive layers can be approximated 
as solid boundaries. To test the reliability of this assump- 
tion, a few cases were run where the initial state consisted 




Fig. 7. Snapshot of a simulation including the 
free interface between two layers. left: solute, 
right: temperature. Pr=1.0, Le=0.01,i?p=1.15, 
Ra* = 6 10^ . See mo vie at j http: / /www.mpa-| 
garching.mpg.de/~henk/doubleJayer.avi 



of two instead of a single step. In some, though not all of 
these runs, the division into two layers remained till the 
end. An example is shown in Fig. [7] for a case with Pr = 1, 
Le = 10-^, Rp = 1.15, Ra* = 6 x 10^ and a resolution of 
500 X 500. Note the approximate (anti-) symmetry of the 
plumes near the interface in the middle, a phenomenon 
known from laboratory experiments. It is caused by the 
continuity of the horizontal velocity across the interface en- 
forced by viscosity. This symmetry becomes less marked in 
simulations at lower Prandtl numbers. 

Fig. |8] shows horizontal and temporal average profile of 
temperature and solute with height. The transition in the 
middle is broad compared with the boundary layers at top 
and bottom. Inspection of the time dependent flow shows 
that this is due to two separate effects. One is the displace- 
ments of the interface by surface waves, which smoothen the 
average gradient without changing the actual fluxes across 
the interface. In addition there is a real mixing effect asso- 
ciated with breaking of the surface waves, but it remains 
localized around the interface. So as to maintain the same 
heat flux, broadening by this mixing has to be compensated 
by a larger amplitude of the transition. The effect is equiva- 
lent to a decrease of the coefficient a in ( 43 1 , and a decrease 
of the effective diffusion rate. 



6. Application to stars 

6.1. The Iow-Pt, low-Le limit 

The astrophysical limit, where the thermal diffusivity is or- 
ders of magnitude larger than the other diffusivities, actu- 
ally simplifies the physics of the double diffusive problem. 
As predicted by relation (381 and in agreement with our 



numerical results, the solute flux is low compared with the 
thermal flux carried by the flow, in this limit. The amount 
of solute which is in the process of being transported from 
the lower to the upper boundary is therefore small at any 
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Fig. 8. Average temperature and solute profiles with height 
z in the double-layer simulation of Fig. [7] 

point in time. The effect of its buoyancy on the dynamics 
of the flow is then only a perturbation. 

To the extent that the solute and viscous boundary lay- 
ers are thin compared with the thermal boundary layer, the 
flow the interface between layers can be approximated as 
a sharp jump. The position of the interface fluctuates due 
to the presence of surface waves, but these waves do not 
transport a net flux (averaged over a wave period, say) . 

The effect of the interfaces is thus equivalent to those 
of solid boundaries, and the flow inside the layer is almost 
the same as thermal convection between plates at the same 
temperature difference. 

Because of this useful circumstance, the heat flux under 
the low-Pr, low-Le astrophysical conditions can be found 
from the equivalent thermal convection problem, and the 
solute flux then follows simply from relation ( 38 ) . For this 
equivalent convection problem, laboratory results and ana- 
lytical estimates are available, and these can then be used to 
extrapolate the numerical results to higher Rayleigh num- 
bers. 



6.2. Extrapolation 

Convection experiments in gaseous Helium at temperatures 
just above the critical point by Castaing et al (1989) yielded 
the fit 



Nut = 0.23 Ra' 



0.282 



(40) 



over the range 10^ < Ra < 10^'*. More recent measure- 
ments by Niemela et al. (2000) using superfluid "^He, over 
the remarkable range 10^ < Ra < 10^^ gave a marginally 
different result: 



Nut = 0.124Ra' 



0.309 



(41) 



This can be compared with the estimate based on a 2- 
dimcnsional argument in S92: 



Nut ~ 0.5 Ral^ 



(42) 



Since the Prandtl number in the laborato ry expe riments 
(of order 0.7) is not far from unity, Ra in (40 41 1 can be 
identified with Ra* in (42). The expressions above are of 
the form Nut = aRa^. As in the above, we subtract 1 



from the Nusselt number, for a marginally better fit at low 
Nut : 



Nut 



1 



^Raf. 



(43) 



To make these things astrophysically useful, we need to 
express Ra* in terms of the superadiabatic gradient (eq. 

Ra* = (V- Va)i^e^ (44) 



where 



R = gH^/K\ 



(45) 



is a function the local thermodynamic state, but not of the 
(still unknown) temperature gradient, and e = d/ H is the 
double diffusive layer thickness d in units of the pressure 
scale height H. 

In the limit e — > 0, this Rayleigh number is equiva- 
lent to that in a Boussinesq calculation. To translate the 
corresponding Nusselt number into an astrophysical flux, 
however, the complication discussed in section [3] has to be 
taken into account. As derived there, the identification (eq. 
17 1 has to be made in terms of the ratio / of convective flux 
Fc to the part i^i-g of the radiative flux that is carried by 
(only) the superadiabatic part of the temperature gradient. 
In terms of the logarithmic gradients, this ratio is: 



F-F, ^ Yl_ 



With ( |43|44|45[ ) this yields 



Nut - 1 = f = (Vr - Va)T^e^ - 1 = aRaf 



) ^ 



(46) 



(47) 



The last equality defines a relation between the layer thick- 
ness d and Ra*, under the assumption of a Nusselt number 
of the form ( 43 ) . With ^ , this then yields the superadia- 
baticity as a function of a. 

The mixing rat e, as fu nction of the layering thickness d, 
follows from eqs. (38 43). Introducing the effective solute 
diffusivity Kge- 

Ksc/«s = Nug = [1 + [KT/Ksy/yRp] (Nut - 1) + 1, (48) 



where (eq. 11 ) 



(49) 



In the limit Nug ^ 1, Nut ^ 1, and assuming the 
second term in the first bracket to dominate, this reduces 
to: 



KSo/'«S = Nug = ( ) ' 



(50) 



This is the same as derived in S92, except that the con- 
tribution of radiation pressure in the equation of state was 
also taken into account there, which yielded: 



KSe = (ktKs) ' I ^ ~ V 



(51) 



where /3 = Pg/ (Pg + Pr) is the ratio of gas pressure to total 
pressure. 
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Fig. 9. Semiconvection for the parameters of a 15 Mq MS 
star. Top left: modified Rayleigh numb er a s a function of 
the double-diffusive layer thickness d (eq. 47), right: Nusseh 



number, bottom left: the corresponding superadiabaticity, 
right: density ratio. Sohd line: using the fit J4l| ), dotted: fit 
(42). The numerical results shown in Fig. jsTlIe just below 
tne solid line in the top right panel, between d/H = 10^'' 
and d/H = 10-3. 



6.3. A 15 Mq star 

We are now in the position to estimate the range of pa- 
rameter values for semiconvection in a star, and which 
part of this range is covered directly by numerical simula- 
tions. Consider the important case of massive stars around 
main sequence turnoff. We use a model kindly provided by 
A. Weiss (model 'fzml5_151'). Characteristic values for the 
physical quantities in the semiconvective zone this model 

0.02, V, 



are g « 10*^ cm/s^. Kg « 1, kt « 3 10*^ cmVs, H Ri 2 10^" 



0.4, Vr - Va 



1. 



Fig |9] shows the 
dependence of Ra*, Nut, V — Va and Rp on e = d/H 
as determined from (471, for two of the power laws above. 
The lower left panel shows how the superadiabaticity in- 
creases with the decrease of convective efficiency at small 
layer thickness. At thicknesses less than 10"* convection 
becomes inefficient at transporting heat, and V saturates 
to its radiative value Vr. The increase of the density ratio 
mirrors this: because of the efficiency of convection at high 
layer thickness and the resulting weak superadiabaticity it 

reaches very high values. 

The effective He-diffusivity from (51) is Kse ~ 10^ 



cm^/s. The mixing time scale over a pressure scale height 
is thus about 10^° 



6.4. Discussion 



As in the case of ordinary convection, there are far fewer 
intrinsic parameters in semiconvection than physical quan- 
tities defining the physical state in a stellar interior. This al- 
lows a significant volume of astrophysically relevant param- 
eter space to be covered by a grid of numerical simulations. 
The low value of the viscosity and constituent diffusivity 
compared with the thermal diffusivity present a limiting 
case that actually simplifies the double diffusive problem 
greatly. In this limit the results become nearly independent 



of the Prandtl number. The remaining parameters can be 
represented by the dimensionless quantities R (eq. 45 1 for 



the thermodynamic state of the plasma, the Lewis number 
Ats/^ Ti the density ratio Rp and the layer thickness d/H 
(eqs. fiTpOl ). 

The sirnulations were all done in 2-D, so tests in 3-D 
will be needed for verification. It is unlikely, however, that 
the results will turn out very different, at least within an 
astrophysical factor of 2. The reason is that in the low-Pr, 
low-Le limit the flow in the layers is almost equivalent to 
ordinary convection between plates. Due to the low solute 
diffusivity, the amount of solute in the bulk of the layer is 
small. It can then be treated as a perturbation, as assumed 
in S92. Known laboratory results for convection at very 
high Rayleigh numbers can then be used to extrapolate 
the numerical results. As shown in section |6] (Fig. |9 1 , this 
makes predictions similar to the simple model used m S92. 
In particular the effective mixing rate is very low. 

Stochastic fluctuations in the flow produce scatter in 
the fluxes of heat and solute, which affect the measured 
averages (cf. Table 1). The accuracy of these averages could 
be improved with longer runs. 

Most of the results presented are based on simulation 
of a single double diffusive layer. As we found in |5.5[ this 
does not capture broadening of the interfaces between lay- 
ers by breaking surface waves (see movie at Fig. [?]). Since 
shallower gradients imply lower fluxes, such a broadening 
process reduces the effective mixing rate. 

The physics determining the thickness of the individual 
double diffusive layers remains uncertain. In the equivalent 
geophysical examples semiconvective zones always consist 
of many very long-lived layers. In the absence of a good 
theory for the layer thickness and its evolution it is not 
possible to translate this flnding to an astrophysical setting, 
however. 

Observations in the east-african rift lakes (e.g. Schmid 
et al. 2010) show that layers first forming at the boundary 
of an expanding double-diffusive zone are always thin, sub- 
sequently growing slowly by a process of merging. This is 
likely to happen in a growing stellar semiconvective zone as 
well. Unclear, however, is how many layers will be left at 
the end of the semiconvective phase. If only a few are left, 
the location of the boundaries between these layers may 
well vary somewhat randomly between stars. This would 
introduce a stochastic element in the internal structure that 
might be related to observed variations in evolution tracks. 



6.5. Conclusions 

The results show that the behavior of double diffusive con- 
vection as observed in geophysics and the laboratory ap- 
plies also to the astrophysical case of low Prandtl and Lewis 
numbers. In particular, the process takes place in the form 
of the characteristic double diffusive layering known and 
theoretically understood since the 1980's. This is in con- 
trast with the models of linear or nonlinear oscillations on 
which most recipes used in stellar evolution are based. The 
results show that the well known 'square root' relation be- 
tween the solute flux and the thermal flux ( 38 1 observed 



in laboratory experiments, holds also in the astrophysical 
case of high thermal diffusivity. This relation is the most 
important factor determining the effective mixing rate. 

By a direct comparison of Boussinesq and fully com- 
pressible simulations we find that the two give equivalent 
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results in the limit of small layer thickness. The Boussinesq 
calculations even capture the essence for layers as thick as 
a scale height. 

We find that in the limit of low Prandtl number the 
results become nearly independent of Pr. This greatly re- 
duces the parameter space to be covered by the grid of 
calculations. 

On the basis of these results wc give fitting formulas 
from which supcradiabaticity and mixing rate can be cal- 
culated as functions of the physical parameters V^, Va, 
R — gH^/n^, Le= ks/kt and layer thickness d. For the 
stellar model used for illustration, however, the effective 
mixing rate is essentially negligible. 
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